Hydrologic Design 4501
Chapter 7: Unit Hydrograpgh and Routing
Ardeshir Ebtehaj
University of Minnesota
1- Linear Reservoir and Convolution
A linear reservoir is a simple and meaningful representation of the watershed response to precipitation events. Continuity equation for a linear reservoir is written as follows:
Figure 1: Schematic of a linear reservoir. A reservoir is linear if we assumed that the storage function
is linearly related to the outflow through the constant K. where
is the outlfow,
is the inflow
(rainfall excess),
is change in the basin storage over time, and K represents a time scale
that captures the residence time of the reservoir. Here, the storage
is the sum the rainfall excess minus the outflow up to time t. Now upon substitution, we obtain and have to solve a first-order ordinary differential equation (ODE):
Multiplying both sides by
, one can get:
Now because of the product rule, we can rewrite the above equation as follows:
,
and thus, we have,
Therefore, we can obtain the outflow of a linear reservior as follows:
1-1 Concept of Convolution
The above solution states that the outflow at any time is a linear sum of the contributions by an initial flow (first term) that goes to zero over time and the integral term (second term) that represents the effects of the reservoir resident time Kon the outflow. The integral term is a convolution integral. Learning the concept of convolution is crucial for understanding the concept of the unit hydrograph.
Mathematical convolution for two generic functions
and
is defined as follows:
Figure 2: Visualization of two examples of the convolution integral. We can see that the convolution of a function
with a kernel
may result in translation (shift) and dilution of
. We will show later on the the rainfall excess can be represented by
and the kernal function or the unit hydrograph is a characteristic of the watershed. The result of the convolution of excess rainfall with the unit hydrograph will the be watershed outflow.
1-2 Instantaneous Unit Hydrograph (IUH)
Now looking back on the integral in the linear reservoir outflow equation, we can see that the inflow
resembles
and the convolution kernal
resembles 
and
. If we define the input rainfall excess as a unit impulse function as follows:
Then, we can define the impulse response function of a linear reservoir as follows:
The function
is often called the instantaneous unit hydrograph (IUH), which is the response of a linear reservoir to a unit excess rainfall with an infinitesimally small duration.
1-3 Unit Hydrograph
Similar to the above formalism for a linear reservoir, a watershed has a response function to a unit pulse of excess precipitation --- albeit a more complex one. This response function is known as the unit hydrograph, which is a time-invariant characteristic of a watershed. More specifically, the unit hydrograph
of a watershed is defined as tthe direct runoff hydrograph (DRH) resulting from a unit pulse of excess rainfall, which is defined as a unit depth of excess rainfall (i.e., 1
, 1
, etc) distributed uniformly over the entire watershed for an effective duration of
. Figure 3: The response of a linear reservoir to a unit impulse (left) resembles the response of a watershed, called unit hydrograph, to a unit pulse of rainfall excess (right), which is a unit amount of rainfall, uniformly distributed over the watershed for duration of
. Note that convolution is a linear operator
which means that we can superimpose the impacts of two pulses of excess rainfall over a watershed through linear combination of their individual outflow responses.
Figure 4: Superimposition of the response functions (i.e., re-scaled and shifted unit hydrographs) to excess rainfall hyetographs. For example, the above schematic shows the response to two pulses of excess rainfall, where
.
Therefore, having the unit hydrograph, one can obtain the outflow of a watershed as convolution of excess rainfall
with the unit hydrograph
as follows:
Below is an example of hydrograph superposition used to obtain the outflow of a watershed in response to a 24-hr storm -- when a 6-hourly unit hydrograph is available.
Figure 5: Example of 6-hr Unit Hydrograph convolution (Credit: The COMET Program).
2- Computation of the Unit Hydrograph
Figure 6: Convolution of excess rainfall hyetograph (ERH) at
intervals with the unit hydrograph of a basin, sampled at
points. The convolution results in a Direct Runoff hydrograph (DRH) at
points, where
.
Since, hydrologic data is not continuous, we must use a discretized version of the convolution equation as follows:
For the example shown in Figure 6, we have
where
,
. Therefore, the above convolution can be expressed as a linear system of equation as follows:
Through convolution, one can obtain the above set of equations through flipping the rainfall excess hyetograph and slide it over the coordinates of the unit hydrograph and then compute the sum over multiplications.
2-1 Forward substitution and least-squares
We can solve the above set of equations through forward substitution. In other words, one can obtain
and then use the second equation to obtain
and so forth. One can arrange all data of the above linear system of equations in a matrix form as follows:
and solve it for
values to obtain the unit hydrograph from the observed direct runoff hydrograph of a precipitation excess event. We can see that in the above linear system of equations we have more equations that unknown. Therefore, calculation of the unit hydrograph is often an overdetermined problem.
The above system of linear equation, which is resulted from a discrete convolution operation can be rewritten in matrix form as follows:
where
is an
column vector,
is an
matrix, and
is an
column vector.
For the provided example, the matrix form is as follows:
The matrix form can be easily implemented into a software library for numerical linear algebra. As previously explained, the equation
is an over-determined linear system, since we have more equations than unknowns. Two of the main methods used to solve this system of equations for
are: 1. Forward/backward substitution: Write out the first or last L equations. Solve the first or last equations and sequentially substitute and solve the remaining equations. In this formalism only L equations are used, which result in multiple solutions to the problem.
2. Least squares: Using linear algebra, we can derive a formula to calculate the least squares estimate of the unit hydrograph ordinates as follows:
where
is the matrix transposition operator, which flips the elements of a matrix over its diagonal elements. Note that by least squares, we mean that the solution has minimum error variance amoung all other possible solutions and the solution is unique. The above matrix formulation of the linear system can be solved efficiently using a technical programing language such as MATLAB. Note that, in the above formulation, we multiply both sides with
to obtain a square matrix
, which might be invertible. We need to note that
is often a rectangular matrix with more rows than columns (
) and cannot be inverted. The matrix
is called pseudo inverse of
. In MATLAB, computation of the unit hydrograph can be simply implelmented using either of the following commands: u = P\q or u = pinv(P)*q. --------------------------------------------------------------------------------------------
Example Problem 7.1: Find the half-hour unit hydrograph using the excess rainfall hyetograph and direct runoff hydrograoh shown in the following table.
Solution: Given
,
and since
thus
.
% Least squares solution for Unit Hydrograph
% Direct runoff hydrograph (DRH)
q = [423 1923 5297 9131 10625 7834 3921 1846 1402 830 313]';
% Excess rainfall Hyetograph (ERH)
p = [1.06 1.93 1.81 0 0 0 0 0 0 0 0];
r = [1.06 0 0 0 0 0 0 0 0];
plot(t',q,'DisplayName','DRH','LineWidth',2)
plot(t(1:length(u)),u,'DisplayName','UH','LineWidth',2)
ylabel('\bf Discharge (cfs)')
bar(t,p,'FaceAlpha',0.3,'DisplayName','Precipitation')
ylabel('\bf Precipitation excess(in)')
--------------------------------------------------------------------------------------------
Example Problem 7.2: Calculate the streamflow hydrograoh for a storm of 6 inch excess rainfall over three hours with hourly rate of p = [2,3,1] using the half hour unit hydrograph computed in the previous example and assume the base below is constant at 500 [cfs] throughout the flood. Check that the total depth of direct runoff is equal to the the total depth of excess precipitation for the watershed area of 7.03 miles square.
Solution:
%computation of the strearmflow hydrograph
u = [404, 1079,2343,2506,1460,453,381,274,173]; % [cfs/in]
Q_direct = conv([0 0 u 0 0],p,'valid'); % [cfs]
Q = Q_direct+Q_base; % [cfs]
V_d = sum(Q_direct)*dt; % [cfs.hr]
V_d = V_d*3600/1; % [ft^3]
Q_direct_1 = conv([0 0 u 0 0],[2 0 0],'valid');
Q_direct_2 = conv([0 0 u 0 0],[0 3 0],'valid');
Q_direct_3 = conv([0 0 u 0 0],[0 0 1],'valid');
plot(t,repmat(Q_base,size(t)),'-ok','DisplayName','Base Flow','LineWidth',2)
plot(t,Q_base+Q_direct_1,'-or','DisplayName','Response to p = [2,0,0] [in]','LineWidth',2);
plot(t,Q_base+Q_direct_1+Q_direct_2,'-og','DisplayName','Response to p = [2,3,0] [in]','LineWidth',2);
plot(t,Q_base+Q_direct_1+Q_direct_2+Q_direct_3,'-ob','DisplayName','Response to p = [2,3,1] [in]','LineWidth',2);
ylabel('\bf Discharge (cfs)')
bar(t,[p, zeros(1,8)],'FaceAlpha',0.3,'DisplayName','Precipitation')
ylabel('\bf Precipitation excess (in)')
2-2 S-Hydrograph
As we mentioned, the UH is a time invariant characteristic of a watershed in response to a unit pulse of excess rainfall with a specific duration
. However, since the UH is obtained through a linear process, we can change a unit hydrograph from one duration to another. To that end, we need to derive the S-hydrograph based on the principle of superposition. By definition, the S-hydrograph is the direct runoff hydrograph of a watershed resulting from a continuous excess rainfall at a constant rate of 1 [cm/hr] or [in/hr] for an infinite period of time.
In other words, the S-hydrograph
is the response of the watershed to summation of an infinite number of unit pulse of excess rainfall with intensity
and duration
-- or a unit rainfall excess over an infinite period of time (Figure 7).
Figure 7: Conceptual representation of the S-hydrograph.
As a result, the unit hydrograph for
duration, with intiger n values, can be obtained as follows:
Figure 8: Derivation of the unit hydrograph with duration
, where
, from an S-hydrograph that is generated from a unit hydrograph with duration
. --------------------------------------------------------------------------------------------
Example Problem 7.3: Use 0.5-hr unit hydrograph in the following table to calculate 1.5-hr unit hydrograph.
Solution:
A shifted version of the hydrogrsph is shown between columns 2 and 3. Column 3 represents the sum of the unit hydrographs multiplied by
[hr]. Column 4 shows the lagged unit hydrograph for
[hr] An excel file solving an example to derive a 1.5-h unit hydrograph from 0.5-hr unit hydrograph is posted in the course website. Column 5 computes the unit hydrograph for
[hr]. % computing 1.5-hr hydrograph from 0.5-hr hydrograph using S-hydrograph
u_0_5 = [404 1079 2343 2506 1460 453 381 274 173 0 0 0];
M = zeros(length(t),length(t));
u_temp = u_0_5(1:end-i+1);
M(i:12,i) = u_temp(1:end)';
s_hydro_lag = [0;0;0;s_hydo(1:end-3)];
u_1_5 = 1/dt_prime*(s_hydo-s_hydro_lag);
2-3 Synthetic Unit Hydrographs
Thus far, the methods we have described only apply for the basin with gage measurements of the streamflow, however, many basins are ungauged. Therefore, engineers have analyzed many unit hydrographs from gauged basins to obtain a universal functional representation for the unit hydrograph that can be generalized to ungauges basins. This representation of the unit hydrograph is often called Synthetic Unit Hydrograph (SUH) that can be applied to ungauged basins. These synthetic methods typically scale a general UH shape or calculate key ordinates based on easily determined characteristics of watershed and excess precipitation.
There are many of these methods, such as Snyder's or Clark's SUH; however, we will cover the very popular NRCS (SCS) dimensionless synthetic hydrograph method. The dimensionless SCS hydrograph and its triangular representation are shown in Figure 9 and tabulated in Table shown as Figure 11.
Figure 9: Dimensionless curvilinear and equivalent triangular SCS unit hydrographs (USDA, 1986)
Input parameters:
: duration of excess rainfall
: time lag between the centroid of rainfall excess hyetograph to the peak discharge
: time to concentration
: time to peak
: recession time
: base time for total duration of the scaled SUH
: peak flow of the scaled SUH.
2-3-1 Concentration time
As we previously discussed the time of concentration (
) for a watershed is the time for a raindrop to travel from the farthest point in the watershed to the outlet. The SCS recommends two methods for calculation of
: - The lag method: defines the time lag as the time in hours from the center of the mass of the rainfall excess to the peak discharge, which can be computed as follows:
where
is the average slope of the watershed in %, L denotes the hydraulic length of the main channel in feet, and
[inches] represents the potential maximum retention that is obtained from the curve numbers (CN). The SCS method suggests the following relationship for the concentration time - The upland or velocity method defines the concentration time as follows:
where
and V are hydraulic length and velocity in [ft] and [ft/s], respectively. The velocity of the upland method can be obtained from Figure 8 based on the average slope of the watershed. Clearly, if we discretize the the stream into k segments, the concentration time is the sum of travel time for different stream segments as we discussed before,
. Figure 10: Velocities for upland method for estimation of
(USDA, 1986). 2-3-2 Time to Peak
Time to peak (
) is the time from the beginning of the rainfall to the time of the peak discharge, which can be expressed as a function of the time lag
as follows:
where
is in
,
is the duration of the rainfall excess in
, and
is the lag time in
. When
is not known, the SCS method suggests and because
, one can obtain the following relationship for approximating the time to peak:
2-3-3 Peak Discharge
From continuity equation, we know that the area under the unit hydrograph is equal to the volume of the excess unit rainfall - example 1 [in] all over the watershed. Here, the unit of
in SCS method is the unit of discharge per depth of rainfall excess such as
or
. Therefore, base on conservation of mass of the depth of the rainfall and the triangular shape of the proposed SCS unit hydrograph, we have
which can be rearranged as
For the SI unit, in the above equation the area should be in
and time is in second to obtain the discharge in
. To convert the above formula to be used for watershed area in square kilometers and hourly time scale, we need to to multiply the above equation by
. Moreover, the SCR method suggests
and thus where
. The
value in English unit is 483.4, where A is in square miles. After obtaining the
and
, the curvilinear coordinates of the SUH can be obtained from the table in the next slide. Figure 11: Ratios for dimensionless SCS synthetic unit hydrograph and mass curve.
2-3-4 Construction of SCS synthetic UH
- Compute the time of concentration (
) using the lag method or the upland-velocity method. - Compute the time to peak
and then the peak discharge
. - Compute the the base time of the SUH as
, which is
for the triangular hydrograph and
for the curvilinear one and the recession time
. - Compute the SUH ordinates and plot them. For the the triangular SUH, the ordinates can be directly computed while for the curvilinear one the dimensionless ratios
shall be used.
Figure 12: The SCS unit hydrographs (USDA, 1986)
---------------------------------------------------------------
Example 7.2: Construct a 10-minute unit hydrograph for a basin with an area of 3.0
and concentration time
[hr].
[min] = 0.166 [hr] (duration of excess rainfall)
[hr] (lag time)
[hr] (time to peak)
(peak flow discharge)Figure 13: The calculated SCS SUH.
4- Flow Routing
Flow routing is the procedure of determining the downstream hydrograph from hydrographs at one or more points upstream.
As water enters a channel, either from hillslopes or upstream channels, the flow velocity shifts the location of the peak flow and surface roughness and storages along the water path attenuate the magnitude of the peak. The flow routing is referred to a class of methods that calculates the redistribution (widening) and translation of upstream hydrograph as a function watershed flow characteristics.
Figure 14: Left: Schematic of redistribution and translation (Chow et al. 1988). Middle: The storage accumulation and release during the flow routing. Right: Conceptualization of channels as storage reservoirs (COMET Program)
There are two main types of routing methods:
1. Lumped (Hydrologic) Routing: Flow is calculated at a particular point as a function of time alone.
2. Distributed (Hydraulic) Routing: Flow is calculated as a function of time and space.
For this course, we focus on the simpler lumped modeling.
Within lumped modeling, there are multiple methods depending on what system we are trying to model. We are going to focus on two common methods for flow routing in reservoirs (level-pool method) and stream channels (Muskingum method).
Derivation of routing methods begins with the continuity equation:
In routing,
is the known inflow hydrograph, while
and
are unknown. Just as in our discussion of linear reservoirs, we need a storage-discharge relationship,
. There are two types of storage-discharge relationships:
- (1) invariable reservoirs: The invariable are those in which there is a one-to-one relationship between the outflow and the reservoir storage. This scenario is valid in short and wide reserviors.
- (2) variable reservoirs: In variable reservoirs, the relationship between the outflow and storage capacity is not unique and follows a loop when the reservoir is filled or emptied. This nonlinear relationship mostly occurs in long and narrow reservoirs where the water velocity is significant and produces backwater effects in the reservoir.
Figure 15: Invariable reservoir (left) that is wide and deep compared to its length. The flow velocity is low and the water surface is almost horizontal. A variable reservoir (right) is a long, narrow and shallow reservoir with high flow velocity such that the water surface profile is markedly curved due to the backwater effect.
Figure 17: Photos of a reservoir, dam and a pond with hydraulic structures as flow controls.
4-1 Level Pool Routing
Level pool routing typically applies to man-made hydraulic structures like reservoirs, dams and stormwater ponds that are common in many watersheds. The key assumption is a horizontal water surface and negligible velocity in the reservoir. This allows an invariable relationship between storage (water surface height) and discharge, because any change in storage results in a uniform change in height of water surface.
Beginning with the continuity equation
, we have
, one can discretize it as follows assuming that the changes in inflow and outflow are linear over a small timestep
: Figure 17: Discretization of the reservoir continuity equation over a small time interval.
The inflow
and
are known at all times while initial storage
and outflow
are specified as initial conditions, so we can re-arrange the equation into knowns and unknowns as follows:
To use this equation, we need the elevation-storage function
and elevation-discharge
relationships of the reservoir to create a storage-outflow function that allows obtaining the outflow at time step
from the left-hand side of the above equation
iteratively.
- The elevation-storage function
is simply constructed from the geometry of the reservoir or pond. - The elevation-discharge
is slightly more difficult to obtain; however, the outflow of storage facilities is often controlled by hydraulic structures such as sluice gates, weirs, pipes, and orifices which typically have well-defined equations based on thr water head (Figure 17).
For example, for an uncontrolled ogee crest spillway, the outflow as a function of total water head on the crest is:
is the effective length of the crest
is the discharge coefficient
is discharge in [cfs].Figure 18: Discharge relationship for vertical-faced ogee crest (US bureau of reclamation 1987). The figure is for English system and
is the water surface upstream from weir drawdown. For a general reservoir, once the storage-outflow function is calculated, as schematically shown in Figure 19, one can calculate the RHS of Eq.(1) and compute the LHS.
Figure 19: Change of storage during a routing period and development of the storage-outflow function for level pool routing on the basis of storage-elevation-outflow curves (Chow et al. 1988).
---------------------------------------------------------------
Example 7.3: Inflow hydrograph is specified at
min. Horizontal area of the reservoir is
acre (43,560
) and
. Given the inflow hydrograph, compute the outflow hydrograph.
Figure 20: The left four columns represent the problem inputs and construction of the storage-outflow curve
. The last six columns represent the inflow and the routed outflow. An excel worksheet is posted in the class website. Solution:
- First we need to construct Q (y-axis, col 2) versus
(x-axis) (col 4) - Claculation of the outflow hydrograph at
:
[cfs] Now we have to use the curve
to obtain
from
:
[cfs] (interpolation) - Calculation of the outflow hydrograph at
:
[cfs]
[cfs] (interpolation) - Computation continues for further time steps ⋯
The whole pocess can be formulated as follows:
- (1) Use the following discretized equation and compute the right hand side form data in the prevoisu time step:
- (2) Conduct a linear interpolation on storage-outflow curve for
for
as follows:
- (3) Given
, compute the right hand side of Eq. 1 for the next time step as follows
- (4) Use the following discretized equation and compute the right hand side form data for the prevoisu time step and continue the iteration ...
Figure 21: The storage-outflow curve (left) and inflow and routed outflow hydrographs for the level-pool routing example.
4-2 Hydrologic River Routing
In order to calculate the routing for streams in a watershed, we can no longer assume that the water surface is level, since there is not negligible velocity in a stream and there are backwater effects. Therefore, we need to define a different storage-outflow equation,
. The Muskingum method assumes that storage in a channel is the sum of a prism and a wedge storage.
- The prism storage is a linear reservoir defined as
where K is a constant representing approximately the travel time through the reach. - The wedge storage is defined as
, where X is a weighting factor ranging from 0 to 0.5 and thus,
This wedge storage represents a flood wave as it propagates through a stream channel as shown in Figure 22.
Figure 22: Schematics of wedge storage and flood wave propagation assumed in Muskingum method (Credit: The COMET Program)
We can discretize this storage equation as follows for time step jth
and thus:
Recall, we previously defined the change in storage as follows:
Combining equations 2 and 3, one can obtain the following flood routing equation:, known as the Muskingum stream routing formulation:
where
. Typically, the streams is broken up into N subreaches when performing this calculation. The U.S. Army Corps of Engineers recommends the following condition to ensure solution stability when we discretize the river reach:
---------------------------------------------------------------
Example 7.4:The inflow hydrograph to a river reach is given below. Determine the outflow hydrograph of this reach, if
[hr],
and
[hr]. The initial outflow is 85
.
Figure 23: The results of the Muskingum routing.